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ABSTRACT 

Context. Mantles of iced water mixed with carbon monoxyde, formaldehyde, and methanol are formed during the so-called prestellar 
core phase. In addition, radicals are also thought to be formed on the grain surfaces, and to react to form complex organic molecules 
later on, during the so-called warm-up phase of the protostellar evolution. 

Aims. We aim to study the formation of the grain mantles during the prestellar core phase and the abundance of formaldehyde, 
methanol, and radicals trapped in them. 

Methods. We have developed a macrosopic statistic multilayer model that follows the formation of grain mantles with time and that 
includes two effects that may increase the number of radicals trapped in the mantles: i) during the mantle formation, only the surface 
layer is chemically active and not the entire bulk, and ii) the porous structure of grains allows the trapping reactive particles. The 
model considers a network of H, O, and CO forming neutral species such as water, CO, formaldehyde, and methanol, plus several 
radicals. We ran a large grid of models to study the impact of the mantle multilayer nature and grain porous structure. In addition, we 
explored how the uncertainty of other key parameters influences the mantle composition. 

Results. Our model predicts relatively high abundances of radicals, especially of HCO and CH3O (1CT 9 - 1CT 7 ). In addition, the 
multilayer approach enables us to follow the chemical differentiation within the grain mantle, showing that the mantles are far from 
being uniform. For example, methanol is mostly present in the outer layers of the mantles, whereas CO and other reactive species are 
trapped in the inner layers. The overall mantle composition depends on the density and age of the prestellar core as well as on some 
microscopic parameters, such as the diffusion energy and the hydrogenation reactions activation energy. Comparison with observations 
allows us to constrain the value of the last two parameters (0.5-0.65 and 1500 K, respectively) and provide some indications on the 
physical conditions during the formation of the ices. 
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1. Introduction 

Among the 150 molecules that have been detected in the inter- 
stellar medium (ISM) so far, a significant number are complex 
organic molecules (here-after COMs), carbon bearing molecules 
with more than five atoms. Several COMs have been observed in 
large quantities for two decades in the w arm and dense hot cores 
of massive protostars dBlake et al.ll 1987b . They have received re- 
newed interest in the last few years after the detection of abun- 
dant COMs in solar-type protostars, spec ifically in hot corinos 
dCazaux et al.ll2003t jBottinelli et al.ll2007l). and in the clouds of 
the Galactic Center dRequena- Torres et al.l l2006). 

Astrochemical models have shown that many COMs can - 
not be produced efficiently in the gas phase. iHorn et all ((2004), 
for example, predicted methyl formate abundances to be less 
than 10~ 10 by considering formation pathways only in the gas 
phase, but this molecule has been detected with abundances up 
to 10~ 6 (see references in the caption of Fig. [TJ. On the other 
hand, several observational and experimental works have high- 
lighted the catalytic behaviour of interstellar grains for chemical 
reactions. First, infrared (IR) observations of protostellar sources 
carried out with space or ground-based telescopes have shown 
that formaldehyde and methanol can be significant components 
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of the grain mantles with fracti onal abundances up to 30% with 
respe ct to water in some cases dGibb et al.ll2004t iBoogert et al.l 
2008). These mantles are believed to be formed during the cold 
and dense prestellar core phase mainly via hydrogenation and 
oxydation of CO and O. Second, several laboratory experiments 
simulating cold cloud condition s have confirmed the effic iency 
of the hydrogenation processes. IWatanabe & Kouchl d2002l) for 
instance have succeeded in producing solid formaldehyde and 
methanol at low temperature (~ 10 K) via CO hydrogenation 
reactions in a CO-H2O ice mixture. 

Modelling the chemistry on the grain surfaces is therefore 
crucial if we aim to understand the formation of COMs on the in- 
terstellar grains. Over th e past 30 years, several n umerical meth- 
ods based on ab initio dAllen & Robinson 1977), Monte Carlo 
dTielens & Hagenlll9 82: Cha rnlev et al.ll 19921) . or rate equations 
dHasegawa et alJl992tlCaselli et al.l!998l) approaches have been 
developed. All these methods give a macroscopic description of 
the problem, meaning that they follow the overall behaviour of 
the mantle rather than the single particle. This method, based 
on the rate equations, though the least precise from a physical 
point of view, is the fastest and allows studies of the chemical 
evolution in objects with varying physical conditions as a func- 
tion of the time. In this context, time-dependent models have 
been developed that predict the abundance of many COMs in hot 
cores/corinos by considering two phases: a dense and cold pre- 
collapse phase during which the grain mantles with simple hy- 
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drogenated species are formed, and a warm-up phase, caused by 
the collapse, during which heavier species, and particularly rad- 
icals, can react toge ther on the grain surfaces before sublimating 
into the gas phase (IGarrod & Her bst 2006; IGarrod etal.ll2008t 
lAikawa et al.ll2008l) . A key point of these models is therefore 
the presence of radicals (OH, HCO, CH3O) trapped in the cold 
mantles. In these models, the radicals are assumed to be mainly 
produced by the U V photodissociation of t he neutral species in 
the grain mantles. IGarrod & Herbstl ([2006) tested this assump- 
tion by introducing reactions between OH, HCO, CH3, CH3O 
to produce methyl formate, dimethyl ether, and formic acid. 
Garr od et alj (|2008) then expanded the chemical network of re- 
actions between radicals to form other COMs such as et hanol, 
glycolaldehyde, o r acetic acid. Finally, Aika wa et al. (2008) and 
Awa d et al] ( l2010l) have included the spatial distribution and evo- 
lution of a collapsing cloud with the aim to estimate the size of 
typical hot corinos such as the well-known source IRAS 16293- 
2422. 

Although the abundances of the simplest COMs, such as 
formaldehyde (H2CO) or methanol (CH3OH), are now quite 
well predicted, astrochemical models still fail to reproduce 
the abundance of more complex ones, such as methyl for- 
mate (HCOOCH3), for example. This molecule is assumed to 
be mainly formed on grains via the reaction between HCO 
and CH3O during the w arm-up phase (lAllen & R obinson 19771: 
IGarrod & Herb st 2006), a reaction in competition with the hy- 
drogenation reaction of CH3O, which leads to methanol. Figure 
Q] shows the observed and the predicted methyl formate to 
methanol gas phase abundance ratios as function of the gas phase 
methanol abundance. The plotted observations refer to the abun- 
dance values estimated in the articles cited in the figure caption, 
where the sizes are either directly estimated by interferometric 
observations or indirectly by other considerations for single-dish 
observations. The plotted model predictions refer to gas phase 
abundances. The curves "Garrod-F" and "Garrod-S" have been 



obtained from Fig. 4 and 6 of IGarrod et alj (2008) respe ctively , 
the curve "AwadlO" comes from Fig. 5 of Awad et alJ ( 12010). 



and the curve "Laasll" from Fig. 5 of lLaas et al. Each 
point of the curves represents a different time or temperature (de- 
pending on the model). In these models, the methanol abundance 
increases with time, except in Awad et al. (2010) who considered 
the destruction of the species in the gas phase. We remark that 
while the methanol abundances may suffer of the uncertainty on 
the source size and H2 column density, the methanol to methyl 
formate abundance ratios are almost unaffected by these uncer- 
tainties. Therefore, unless the observed methanol abundance is 
always lower than 10~ 7 , which is certainly not the case, no model 
can reproduce the totality of observations. Specifically, the ob- 
served ratios are roughly independent of the methanol abun- 
dance and are about 0. 1 and 1 in hot cores and hot corinos re- 
spectively, regardless of the telescope used for the observations. 
State-of-art astrochemical models consequently underestimate 
the methyl formate to methanol abundance ratios and always 
predict a decreasing ratio with increasing methanol abundance 
(because methyl formate is in competition with methanol for- 
mation). 

Motivated by this unsatisfactory situation, we developed a 
new model, called GRAINOBLE, whose ultimate goal is to sim- 
ulate the synthesis of COMs on the grain surfaces. In this article, 
we focus on the first phase, namely the formation of the mantle 
during the dense and cold phase of the pre-collapse, the prestel- 
lar core phase. Our model differs from most previous published 
models because it treats the multilayer and porous structure of 
the iced mantles. In the mentioned previous models (Garrod & 



x 

o 

o 

X 

o 

X 



10 l 



10- 



X 

o 
o 
o 

g 10- 4 

X 



10"' 



X 

. 

v x 


1 

X 
O 

\**x 


1 

X X 

x X 


1 

Hot corinos 
Hot cores 

+ 1 s yr 

+7.1 s yr 

X X 
X 


e^r 

. Laasl I^^Cv 


Awad10\ 




\^Garrod08-F 




' 


^— > 

Garrod08-S\ 
' 





10- 



10' 



1 Q' 7 1 f 
X(CH 3 OH) 



10" : 



1(T 



Fig. 1: Gas phase methyl formate to methanol abundance ra- 
tio as function of the gas phase abundance of methanol. 
Observations of hot corinos with single-dish telescopes are 
repres ented by red cross es: NGC1333 - 4A by iBottinelli et al] 
J2004h . IR AS 16293 bv ICazaux et all J20031 NGC1333-2A 
and -4B by IBottinelli et alJ (120071) . Observations of hot cores 
with single-dish telescopes are repre s ented by blue crosses: 
G34.3+0.15 bv iMehringer & Snvderi (fl996h . Sgr B2(N) and 
SgrB2 (M) bv lNummelin et all (12000170 327.3-0.6 bv lGibb et all 
(2000), OMC1 by Sutt on etalJ (11991 . G 34.3+0.2, Sg r B2(N) , 
DR23(OH), W51, Orion Hot-Core by llkeda et alj (I200TI) . 
AFGL25 91, G24.78, G75 78, NGC6334, NGC7538, W3(H20), 
W33A by Bisschop et alj d2007l) . Observations of hot cores with 
int erferometers are r e presen ted by blue diam onds: G34 . 3+0.1 5 
by iMacDonald etal] (119961). G19.61 - 0.23 b y lOin etalJ d2010b . 
Or ion KL, G29.9 6 by Beuth er et al.1 j2009h . and G47.47+0.05 
bv iRemiian et al.1 (|2004). The red curves report the predictions 
of the methyl formate to methanol abundance ratio appropriate 
to the hot corinos case by Garrod et al. (2008, their Fig. 6) and 
Awad et al. (2010, their Fig. 5). The blue curves report the pre- 
dictions appropriate to the hot cores case by Garrod et al. (2008, 
their Fig. 4) and Laas et al. (201 1, their Fig. 5). The arrows rep- 
resent the direction of the time in each model and the ticks refer 
to timescales. 



Herbst 2006; Garrod et al. 2008; Aikawa et al. 2008; Awad et al. 
2010; Laas et al. 201 1), COMs are synthesized during the warm- 
up phase thanks to the increased mobility of radicals, formed on 
the grain surfaces during the previous cold phase. The abundance 
of radicals in the mantles during the cold phase is, therefore, a 
key point. One possibility, explored by the above models, is that 
the radicals are formed by the UV photodissociation of frozen 
species. However, since the impact of the UV photons on the 
surface chemistry is poorly understood (see the detailed discus- 
sion in §2.1), we explore here the possibility that radicals are 
synthesized even without photolysis and are trapped in the grain 
mantles, thanks to the multilayer mantle treatment. 

The article is structured as follows. In the next section, we 
review the models available in the literature, outlining their lim- 
its and pointing out possible improvement. The concluding para- 
graph of the section describes the improvements proposed by our 
new model GRAINOBLE. In section [3] we describe the techni- 
cal details of GRAINOBLE. Section|4]reports the results of the 
new computations, discussing a reference set of parameters and 
a large grid of runs where several input parameters are varied. In 
section [5] we compare the obtained results with previous models 



2 



V. Taquet et al.: The GRAINOBLE model 



with the twofold goal to validate our code and discuss the differ- 
ent results from the different adopted physical models and con- 
ditions. Section [6] compares the GRAINOBLE predictions with 
the observations of ices and discusses the constraints on the var- 
ious parameters of the model. Finally, §|7]concludes the article, 
emphasising the main results and perspectives. 

2. Limits of the existing models and need for 
improvement 

2.1. The physics of grain surface chemistry 

a ) Multilayer vers bulk chemistry 

One of the main limits of many published models is that they 
do not distinguish the chemical processes that occur in the man- 
tle bulk and on the surface. In these models, species that are 
buried in the mantle can keep diffusing in the mantle and reacting 
with others even if they are covered by new layers. Furthermore, 
particles landing on the surface on top of the grain mantle can re- 
act with all particles, regardless of their depth within the mantle. 
Thus, heavy reactive particles such as radicals can continue to 
react with landing hydrogen atoms even though the radicals are 
buried deep in the mantle bulk. The "bulk chemistry method" 
may therefore (grossly) underestimate the abundance of radicals 
in grain mantles. 

In contrast, the difference in the chemical behaviour of the 
mantle and of the surface may lead to the formation of non- 
homogeneou s grain ices as suggested by the IR observations. 
iTielens et alJ {l991), for instance, by observing two features in 
the CO bands, suggested a polar mixture enriched by water in the 
bottom of the mantle covered by a non-polar one mainly com- 
posed of CO located on the outermost layers. In addition, ex- 
periments have shown that UV and cosmic rays irradiation can 
alter the chemistry of astrophysically relevant ices composed of 
water, methanol, and ot her CO-bearing molecules even in dark 
cold c loud conditions (Gera kines et al.l[l996t iBennett & Kaiserl 
|2007|) . However, unlike highly energetic cosmic rays, which can 
cross through the entire grain mantle, UV photons only pen- 
etrate a limited number of monolayers, depending on the ab- 
sorption cross-section of molecules that constitute the ice (see 
iGerakines et al.ll2~00ll for an illustrative discussion). Modelling 
the photolytic processes in ices requires, therefore, a distinction 
between the surface and the bulk of the mantle. 

The desorption of grain mantles during the warm-up phase 
depends on their chemical composition. By carryi ng out temper- 
ature programmed desorption (TPD) experiments. ICollings et al.l 
(2004) demonstrated that CO ices desorb at ~ 20 K, whereas 
methanol and water ices desorb at temperatures higher than 
100 K, for instance. Numerical values of desorption energies 
were deduced from these experiments and were incorporated in 
macroscopic rate eq uation models. Consequently, the model of 
iGarrod et alJ (|2008) predicts a total desorption of the CO reser- 
voir at 25 K, whereas more refractory species are desorbed at 
temperatures higher than 100 K . Actually, a significant amount 
of volatile ices can be trapped with highly polar molecules, such 
as water, in the mantle bulk. Thus, the desorption of such species 
is limited by the gradual evaporation of the mantle main matrix, 
which depends on other molecules. 

The continuous time rando m walk (CT R W) M onte- 
Carlo method, i ntrod uced by Chang et al. (2005) and 
Cupp en & Herbstl (120051) . strongly differs from the macro- 
scopic methods described above. This "microscopic" model 
considers the actual hopping of a discrete number of particles 
from one site to another but, unlike previous Monte Carlo 



methods, it takes into account the actual position of particles 
so th at the spatial distr i bution is al so included in the si mula- 
tion. ICuppen & Herbstl d2007l) and ICuppen et all d2009l) have 
explained the growth of ice monolayers in dense clouds and 
showed a differentiation of the chemical com position in the 
mantle. More particularly, Cuppe n et al.l (12009) showed that un- 
der dense cloud conditions, CO can be trapped in the mantle by 
the accretion of new species above it before reacting with other 
particles. They also highlighted the survivability of very reactive 
particles in the mantle bulk by this process. Unfortunately, this 
method can only study simple chemical systems because of the 
CPU time required by the Monte-Carlo algorithm, and a direct 
lin k with the gas phase i s not possible at present. 

lHasegawa & Herbstl (Q993b) have attempted to treat sepa- 
rately the chemical processes within the bulk and on the surface 
with a three-phase model. They considered the formation 
of an inert mantle caused by the accretion of particles onto 
it. They found that the differences between the three-phase 
and the two-phase models are small for stable species but 
become non-negligible for reactive species at long (> 10 s yr) 
timescales. Furthermore, the values of desorption and diffusion 
energies used in 1993 are lower than the values measured by the 
experiments carried out during the past decade. Because higher 
desorption and diffusion energies values tend to decrease the 
rate of the chemical processes on grain surfaces, the impact of 
the three-phase model on the frozen radical abundances would 
be significantly higher with the new values. However, at the 
time of the submission of this article, this three-phase model has 
not been pursued in other works. 

b) Porous grain surfaces 

Astrochemical models used for predicting the COM forma- 
tion have so far only considere d perfectly sp herical and smooth 
grains. However, observations (Mathis 1996), theoretical studies 
(lOssenkopl 11993b lOrmel et alJ 120091) . and analyses of solar 
system bodies have provided evidence for the fl uffy and porous 
structure of interstellar grains. For example, Mathis (1996) 
needed to use interstellar grains with a substantial fraction 
of vacuum to fit the observed extinction curves. Numerical 
simulations of dust coagulations carried out by Ossenkopf 
(1993) from a MRN-like grain size distribution have shown an 
increase of the size and the porosity of the dust in molecular 
cloud conditions. This fluffy structure of interstellar grains may 
alter the surface chemistry because grains can have micro-pores 
of small apertures compared to their volume in which species 
can be trapped. Thus, porous grains can increase the reactivity 
of chemical reactions by trapping volatile particles such as 
atomic hydrogen in their pores. 

c) Photolysis 

As described in the introduction, the new generation of 
COMs formation models predicts that COMs are synthesized 
during the warm-up phase thanks to the increased mobility of 
radical s, formed on the grain surfaces during the previous cold 
phase dGarrod & Herbsdl2006tlGarrod et al.ll2008UAikawa et al.l 
120081 lAwad et al.l2010tlLaas et al.l201 ll) . The abundance of rad- 
icals in the mantles during the cold phase is, therefore, a key 
point. In these works, UV photons from cosmic rays play a ma- 
jor role in synthesising radicals during this phase. However, as 
for any model, a series of inevitable assumptions are adopted to 
treat the process. 

The first important assumption regards the photodissociation 
rates on the ices. They are assumed to be t he same as those in 
the gas phase, computed by Sternberg e t al.l (ll987l) : lGredel et all 
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(1989). However, there are various reasons to think that this may 
be a serious overestimate. First, grain surfaces could absorb part 
of the UV photon energy, ending up with a lower dissociation 
rate. Second, once the molecule is broken, if it is broken, being 
not in the gas phase but on a surface with almost no mobility, the 
photoproducts may recombine almost instantaneously. In addi- 
tion, some product s can even have eno ugh energy to sublimate 
(see for example I Andersson et al.ll2006l) . Thus, the products of 
the photodissociation are likely different from those in the gas 
phase, and the overall formation rate of the assumed products 
(and radicals) may be severely lower than assumed. 

A second point is related to the flux of CR induced UV pho- 
tons in dense clouds. The exact value is fairly uncertain, be- 
cause it largely depends on the primary cosmic ray energy spec- 
trum, and on the grain exti nction cross section, which are also 
uncertain (see for example iPadovani e t al. 2009). For example, 
IShen et aTl d2004l) have shown that the typical uncertainty on the 
low-energy cosmic ray spectrum leads to a variation of the UV 
flux of more than one order of magnitude. 

Another important assumption regards the penetration of the 
UV photons through the bulk of the mantle. Unfortunately, there 
exists almost no observational or theoretical ground for that, and 
it is not obvious ho w realistic this assumpt ion is. Indeed, the ex- 
perimental work by Gera kines et alj ((2000) suggests that the UV 
photons can only penetrate a limited number of layers, depend- 
ing on the optical properties of the ice. 

Finally, even the branching ratios of species caused by the 
photodissociatio n in the gas phase are very poorly known (e.g, 
lLaas et al]|201 lb . In this cont ext, various authors omit the pho- 
tolysi s in their model (e.g, IChang et all l2007t ICuppen et alJ 
120091) . 

2.2. The computational methods 

The rate equations method allows us to study the evolution of 
the chemical composition of grain mantles via one temporal dif- 
ferential equation for each species. This approach allows us the 
use of complex chemical networks that include photolytic pro- 
cesses involving hundreds of species, and allows us to directly 
link the gas phase and the grain surfaces. However, rate equa- 
tions are based on the evolution of densities of species in gas- 
phase and on grain surfaces, which are average values. Therefore 
this approach can be inaccurate when the number of particles on 
the grain becomes low because of the finite surface of grains. 
ICaselli et al.l £1998) took into account the discrete nature of the 
individual composition of grain mantles by mod ifiying the reac- 
tion rate coefficients, while iGarrod et al] (120081) we nt further by 
modifying t he functiona l form of the reaction rates. Biham et al.l 
d2001b and lGreen et ail f 2001) introduced a stochastic method 
based on the resolution of the master equations (ME), in which 
the system solves the time derivatives of the probability of each 
discrete mantle composition. In spite of its good accuracy for 
grain surfaces, this method can only be used for small chem- 
ical networks since the number of equations grows exp onen- 
tially with the number of species . ILipshtat & Bihaml (12003) have 
thus introduced another method based on the time derivatives 
of the moments to study the f ormation of molecular hydrogen, 
while iBarzel & Bihaml d2007l) have extended this method for 
more complex chemical networks. The comparisons between the 
methods introduced here show that the rate equations method 
tends to overestimate the reaction rates and the abundance of the 
reactants, especially when they are in low quantities. 



2.3. Need for improvement: the GRAINOBLE model 

Although the past few years have seen a huge improvement 
in the modelling of grain surface chemistry, models are not 
yet able to accurately reproduce observations (e.g. Fig. [TJ. 
Therefore, something is still missing in those models. In the 
following, we review what are, in our opinion, the areas where 
improvements are needed and possible. Our new grain surface 
model, GRAINOBLE, has been developed to include the 
improvements described below. 

a) Multilayer structure of the ices 

The first obvious limitation of the Herbst, Garrod, and 
collaborators class of models, the present state of the art, is 
that they do not take into account the multilayering formation 
of the ices, whic h leads to the diffe r entiation observed in the 
interstellar ices dTielens et alj I199U iPontoppidan et alj 120031) 
that is also predicted by the microscopic Monte-Carlo (MC) 
models of Cuppen and collaborators, and to the potential 
survival of reactive particles in the mantle bulk. On the other 
hand, the microscopic MC treatment is too cumbersome in 
computer time and difficult, if not impossible, to apply to model 
realistic cases. GRAINOBLE uses the rate equation method 
developed by Herbst and collaborators (which is at the base 
of the Garrod models) in a way that permits us to follow the 
multilayering structure of the ice. Therefore, it benefits from the 
advantages of the rate equation method (computer speed) and of 
the microscopic MC (multilayer) approach. 

b) Porosity of the grains 

Laboratory experiments, numerical simulations of grain 
coagulation, mantle formation, and astronomical observation s 
all demonstrate that interstellar grains are porous (Jones 2011). 
As previously discussed, the presence of pores creates traps for 
the molecules and atoms on the grain surfaces, leading to an 
enhancement of reactivity on the g rains. GRAINOBLE inc ludes 
pores, following the treatment of iPerets & Bihaml (HOO©- The 
number and area of the pores are treated as parameters (section 
®. 

c) Multiparameter approach 

Many parameters in the microphysics of the problem have 
highly uncertain values, like the various activation barriers of 
key chemical reactions, or the barriers against diffusion of ad- 
sorbed particules. Additionally, several macroscopic parameters 
are also either uncertain or vary depending on the astronomical 
object. In particular, it is worth mentioning that the densities and 
the temperatures of the prestellar cores depend on their age and 
their surroundings. 

3. The GRAINOBLE model 

3. 1 . General description 

Following the work of lHasegawa et alj (ll992l) . four main pro- 
cesses occuring on grains are taken into account. 

1) Gas-phase particles can accrete onto the grains, which are 
considered spherical. For a given gas-phase species i, the ac- 
cretion rate is a function of the thermal velocity of the gas- 
phase species v(i), of the cross section of the grains cr{ad), of 
the density of the grains «</(a</) where a^ is the grain diame- 
ter, and of the sticking coeffici ent S (i) . We as sumed the stick- 
ing coefficient estimated by Tielens (2005) for atomic hy- 
drogen and a sticking coefficient equal to 1 for heavier parti- 
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cles, b ased on the TPD experimental work of Bisscho p et alJ 
(2006). For dark cloud conditions, particles constituting the 
ice are bound mainly via the physisorbed van de Waals inter- 
actions. 

2) Once the particles are stuck onto the grains, they can dif- 
fuse along the surface via thermal hopping according to 
the Boltzmann law, which give s the hopping ra te Ri,„ p . 
Following experimental work of iKatz et al.l d 1999b . we ne- 
glected tunelling diffusion. 

3) Physisorbed p articles can react o nly via the Langmuir- 
Hinshelwood (Hinshelwood 1940) process, in which two 
species can react when they meet in the same site. The reac- 
tion rate R, is given by the product between the diffusion rate 
(i.e, the number of times per second that a species sweeps 
over a number of sites equal to the number of sites of the 
layer) and the probability of the chemical reaction, which is 
a function of its activation energy. 

4) Surface species can desorb only via thermal processes, i) 
the "classic" thermal process caused by the thermal bal- 
ance of the grain, if) the cosmic ray induced desorption pro- 
cess in whi c h cosm ic rays heat the grains, as suggested by 
Lege r~et~aT1 (fl 985). whose heating rate was computed by 
Hasegawa & Herbstl d!993al) . The sum of these two thermal 
desorption rates gives the total evaporation rate R ev . 

A full list of the symbols used in this work is provided in the 
appendix. 

3.2. Porosity 

To model the impact of the porous structures of interstellar 
grains, we introduce t wo types of sites: the no n-porous, and the 
porous sites, following lPerets & B iham (2006). The non-porous 




Fig. 2: Schematic view of a portion of the mantle layer as mod- 
elled by GRAINOBLE. All pores are assumed to be square, to 
have the same size, and no walls (see text). The edges sites, that 
connect the non-porous surface and the pores are depicted in 
blue. In this example, each pore is constituted by 16 sites, of 
which the 12 blue sites are edge sites, and the fraction of area 
occupied by the pores is 0. 1 . 

sites form a smooth surface that is perfectly regular, where par- 
ticles interact directly with the gas phase. The gas phase species 
can accrete onto the grain, and the adsorbates can desorb into the 
gas phase. The porous sites correspond to sites where no direct 
accretion and desorption are possible. They are filled up only 
by diffusion from non-porous sites at the edge of the pores. The 
porous sites cover a fraction F por of the grain surface. For the 



sake of simplicity, all the pores are assumed to be square and to 
have the same size, constituted by N pore sites, where N pore > 4. 
The fraction of the grain surface occupied by the pores, F por , is 
poorly known because, to date, no experimental measurements 
or numerical estimates have been attempted. Therefore we ran 
models with four values for F por : 0, 0.3, 0.6, and 0.9, to cover a 
large range of possibilities, 

The fraction of the edge porous sites, F e d, on the grain de- 
pends then on F por and N pore as follows: 

yj N pore - 1 

F ed = 4-F por ^JL . (i) 

™ pore 

A posteriori, the actual value of N pore does not influence the re- 
sults significantly because what matters is mostly the fraction of 
area occupied by the pores and not the size or the number of 
pores. In the following, we assumed N pore = 9. 

The approach that we used for modelling the grain porosity 
is sketched in Fig. [2] 

3.3. Multilayer approach 

In order to distinguish the chemistry on the surface and in the 
bulk of the mantle, we used a multilayer approach. In practice, 
the code follows the chemical processes on a single layer, the one 
at the surface of the mantle. In this layer, particles can accrete, 
diffuse, react with each other and desorb, depending on whether 
they are on non-porous or porous sites. The mantle layer growth 
is treated at the same time as the chemistry. Specifically, at each 
time t, the chemistry is followed by solving the set of differen- 
tial equations described in §3.9 between t and t + At and a new 
number of particles on the grain is computed alt + At. The layer 
is considered chemically inert as soon as the number of particles 
that are on the layer is equal to the number N s (cid) of (porous 
plus not-porous) sites of the layer, where N s {ad) is given by the 
ratio between the grain surface area and the site area d 2 s (see 
§I3.71 >. Because the porous sites are populated only via diffusion 
(§3.2), they fill up more slowly than the non-porous sites. It is 
therefore possible that the layer is considered inert when some 
pores are still empty. In this case, the number of non-porous 
sites is artificially increased. Although this is somewhat arbi- 
trary, a posteriori this assumption has almost no impact on the 
results, as described in §4.3. When the layer becomes inert, the 
code memorizes its composition and a new and reactive layer is 
started. Note that no chemical exchange is allowed between lay- 
ers. Finally, the code also takes into account the growth of the 
grain size and of the number of sites of the new layer by assum- 
ing that the thickness of the layer is equal to the site size d s . All 
layers have, therefore, the same thickness, equal to the site size. 
Figure |3]gives a schematic view of this approach. 

3.4. Chemical network 

We ad opte d the network based on the work bv lAllen & Robinson! 
(119771) and lTielens&Hagenl (Il982l) . updated for the new theo- 
retical and experimental results and adapted to model the for- 
mation of the main constituents of intestellar ices. As previously 
mentioned in the introduction, observations of ices have shown 
that grain mantles are mainly composed of O- and CO- bear- 
ing molecules such as water, CO, CO2, or methanol. Molecules 
such as NH3 and CH4 produced by the hydrogenatio n of N and 
C can also be found, but only in smaller quantities (Gib b et all 
2004). In this work, we therefore neglected the formation of 
ammonia and methane. In addition, for the reasons explained 
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Diffusion 



L-H reaction 



layer 



Fig. 3: Schematic view of the multilayer approach adopted by the 
GRAINOBLE code. The inner layers of the mantle are chemi- 
cally inert, whereas the layer at the surface is chemically reactive 
(see text). 



in §2.1, we also do not consider photolytic processes. Our 
model, therefore, considers the formation of stable species (e.g, 
formaldehyde and methanol) and radicals (e.g, HCO, and CH3O) 
via hydrogenation of frozen CO. Note that we did not con- 
sider CH2OH because, in absence of photolytic processes, it 
is formed via the hydrogenation of H2CO. Theoretical studies 
have shown that the reaction H2CO + H — > CH2OH has an ac- 
tivation barrier (10.5 kcal/mol) more th an twice the activation 
barrier leading to CH 3 Q (4-6 kcal/mol) dSosa & Schlege]||1986t 
IWalcrJl993l:IWoonl2002h . We therefore neglected CH 2 OH in the 
present study. For the H2O formation, we adopted t he new reac- 
tion scheme based on the e xperi mental works by iDulieu et alj 
d2010l) . Hoppolo et alJ d2010h and lCuppen et alJ (l2010h . who eT 
timated the energy barriers of several reactions. The recombi- 
nation reaction of hydrogen atoms, forming molecular hydro- 
gen, is also taken into account. Since H? has a low binding 
energy with respe ct to water ice ([Tielens & Ailamandolal ll987; 
iPeretset al.ll2005h . and since its reaction of formation releases 
chemical energy dHollenbach & Salpeteri [l970h . we assumed 
that the H2 molecules are desorbed in the gas phase as soon 
as they are formed. Finally, IR observations have shown that 
CO2 may be a major component of ices with abun dances reach- 
ing 30 % with respect to to w ater in some cases (Dartoi s et alJ 
1999; Pontop pidan et al.ll2008l) . However, the formation path- 
ways of carbon dioxide, either by a direct surface route or by 
UV photolysis, in such cold conditions are still not well under- 
stood. Accordingly, we decided to neglect the formation of this 
molecule. 

Table[T]lists the chemical species considered in our network 
with their binding energies relative to an amorphous water ice 
substrate. Indeed, interstellar grains located in prestellar-cores 
are already covered by a polar ice mainly composed of water, 
formed during an earlier phase when atomic abundances are still 
high, before the depletion of CO. The network used for the grain 
surface chemistry is reported in Table[2] 

3.5. Gas phase initial abundances 

The observed and predicted gas-phase abundance of atomic car- 
bon is much lower than CO or atomic oxygen abundances. 



Table 1 : List of species and binding energies assuming that the 
ice matrix is mainly water ice. 



Species 


E b (water ice) (K) 


H 


450" 


O 


800* 


CO 


1150 c 


OH 


2820 d 


H 2 


5640' 


o 2 


1000 c 


O3 


1800 f 


H0 2 


5570 c 


H 2 2 


6300 / 


HCO 


1600 f 


H 2 CO 


2050 f 


CTLO 


5080 ' 


CH 3 OH 


5530 



Refer ences. (a) iHollenbach & SalpeteJ dl970h: w iTielens & Hagenl 
(fl982l); (c) ICollings et all j2004h; w [Avgul & Kiselevl dl970h : (e) 
ISpeedv et all dl996h : ^' iGarrod & Herbstl d2006h 



Table 2: Surface reactions and activation energies. 



Number 


Reaction 


E a (K) 


1 


O + O — > 2 





2 


+ o 2 — > 6 3 





3 


H + H — > H 2 gas 





4 


H + O — > OH 





5 


H + OH — > H 2 





6 


H + o 2 — ^ H0 2 





7 


H + H0 2 — ^ H 2 2 


(36%)" 


8 


H + H0 2 — ^ OH + OH 


(57%)" 


9 


H + H 2 2 — > H 2 + OH 


1400* 


10 


OH + OH — > H 2 2 





11 


H + O3 — > 2 + OH 





12 


H + CO — > HCO 


400 c <E a < 2500 rf 


13 


H + HCO — > H 2 CO 





14 


H + H 2 CO — CH 3 


400 r <E a < 2500 rf 


15 


H + CH 3 — CH3OH 






References. < n) bran chi ng ratios measured by Cuooen 
Kle mm et al] dl975h : (c) lFuchs e t al. (2009); (d) Woon (200! 



et al J 
1002) 



d2010t) : 



Therefore, we considered only the accretion of gaseous H, O, 
and CO onto the grain mantles. 

In dark cloud conditions, the abundance of atomic hydrogen 
in the gas phase results from a balance between the cosmic rays 
ionization of H2 on the one hand and the accretion on dust grains 
(which react to form H2 or heavier iced species) on the other 
hand. At steady state, the density of H is therefo re given by the 
ratio between these two processes dTielensll2 005i. as follows: 



»(//) = 



2.3 ■ Ccr ■ n(H 2 ) 
2 ■ v(H) -cr d -X d - n(H 2 ) 



(2) 



where v(H) is thermal velocity of atomic hydrogen. The result- 
ing H density is independent of the density and equal to 1.2 
cm -3 at 10 K, for a cosmic rays ionization rate ( of 3 ■ 10~ 17 
dCaselli et al.lll998l) . a grain abundance X d of 1.33 • 10 12 rela- 
tive to H nuclei and for a grain radius r d of 0. 1 /mi (<Td = n • r 2 ,), 
giving n d ■ (Jdlnn — 4.2 • 10~ 22 cm 2 . Observations and astro- 
chemical models show that CO is a very stable species and the 
most abund ant molecule af t er H2 in dark clouds. Following the 
estimates of Frerki ng" et alJ (Qj?82), we assumed an initial abun- 
dance of CO of 4.75 • 10~ 5 with respect to H nuclei. This abun- 
dance decreases with time because of the freeze-out onto the 
grain mantles. Finally, the abundance of atomic oxygen in the 
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gas phase remains poorl y constrain e d. Ob s ervations of the [Oi l 
63 um line obta ined by ICaux et al l (Il999t) : IVastel et all (|2000); 
iLis et al.l d2001h with the Infrared Space Observatory (ISO) have 
measured abundances higher than 10" 4 , suggesting that all the 
gaseous oxygen is in atomic form. But astrochemical models 
predict an atomic oxygen abundances of 6 • 10" 5 or less. The 
initial abundance of O is therefore considered as a free parame- 
ter and can vary between 2 • 10" 4 and 2 ■ 10" 5 . 

3.6. Physical conditions 

As supported by observational and theoretical arguments, the 
physical properties of prestellar core s vary substantially with 
the age and the mass o f the object (see ldi Francesco et al.ll2007t 
Bergin & Tafalla 2007). The densities n H evolve typically from 
~ 10 4 to ~ 10 6 cm~ 3 , while the gas kinetic temperature (equal 
to the grain temperature in most cases) varies between 7 and 20 
K, depending on the type of the object. Indeed, massive stars 
are suspected to be formed from warmer prestellar cores than 
low mass ones. Consequently, we included several values of the 
density and the temperature, listed in Table [3] in our parameter 
grid. 

A dditionally, as show n by numerical simulations (Ossenkopf 

19931; lOrmel et al.ll2009l) and observations (IStepnik etalj |2003; 
Steinacker et al.l 1201 Ol) . the coagulation of grains occuring in 
molecular clouds and in prestellar cores tend to increase their 
size. That is why the grain size distribution, which is assumed 
to follow a broad n early power-law in diffuse interstellar clouds 
(Mat his et al.ll 19771) . tends to evolve with time, the bigger grains 
dominating most of the dust mass. We assumed all interstellar 
grains to have the same diameter a,;, and we considered three 
values: 0.1, 0.2, and 0.3 /mi respectively. The grain abundance 
Xj changes in accordance with the grain size, to keep a dust-to- 
gas mass ratio of 1 %. 

3. 7. Other key parameters of GRAINOBLE 

In addition to the above described parameters, astrochemical 
models of grain chemistry are based on several other parameters 
related to the particle-surface, surface particle-particle interac- 
tion and surface properties. We list below three additional key 
parameters of the GRAINOBLE model whose values are poorly 
constrained because of the difficulties involved in measuring or 
computing them. We considered them as free parameters in a 
range defined by the different values found in the literature. 

a) Diffusion energy to binding energy ratio Ed/Eb 

Theoretical calculations of barriers against diffusion for 
physical adsor ption Ed on perfec tl y sm ooth surfaces have been 
carried out by iJavc ock & Parfitt (1986) and have shown that 
Ed is typically about 30 % of the binding energy Ej,. However, 
"real" surfaces, which sh ow defects, irregularitie s, or steps, tend 
to increase this ratio (lEhrlich & Huddal [1966). Accordingly, 
several experimental studies have been carried out to constrain 
the diffusion barrier on astrophysically relevant surfaces. 
By fitting experimental formation of molecular hydrogen on 
refractory bare surfa ces and on amorp hou s water ices with a 
rate equation model, iKatz et al.l ([1999) and iPerets et all d2005l) 
experimentally found hi gh Ed/Eb rati o s (0-7 7 and ~ 0.85 
respec tively ). However, IColl ings et al] (120031) . lUlbricht et al.l 
(2002j), and iMatar et all (120081) estimated lower energy ratios 
(of about 0.5) for CO and D on water ice, and for atomic oxygen 
on carbon nanotubes. Clearly, the energy Ed/Eb ratio is highly 



uncertain and strongly depends on the composition and the 
structure of the substrate. Three values of the Ed/Eb energy 
ratio have been considered here: 0.5, 0.65, and 0.8. 

b) Activation energy E a 

The activation energy values of the CO and H2CO hy- 
drogenation reactions are also very uncertain. Theoretical 
studies of these reactions on an ic y mant l es via calc ulations of 
quantu m chemistry carried out by IWoonl (|2002) and Goumans 
(120111) . showed that the activation energies are about 2000-2500 
K, and depend slightly on the presence of water molecules 
around the reactants. In contrast to these theoretical values, 
experime nts show that these reaction s take place a t low tem- 
perature (IWatanabe & Kouchil |2002|) . iFuchs et all d2009h . by 
modelling the experiments of a hydrogen deposition on a CO 
ice with a CTRW method, found values of about 400-500 K 
for both reactions. We therefore consider these barriers as a 
free parameter. Three values of activation energies, assumed 
to be equal for the H+CO and H+H2CO reactions, have been 
considered: E a = 400, 1450, 2500 K. 

c) Site size d s 

The surface density of sites s, i.e. the number of sites per 
surface unit given in sites x cm" 2 , depends on the structure of 
the surface and its composition. Most astrochemical models 
assume a value of s = 10 15 cm" 2 which has been measured 
experimentally by iJenniskens et al.1 ([1995) for a high-density 
amorphous water ice. By assuming that the distance between 
two sites d s is constant along the surface, d s can be easily 
deduced from the surface density s (d s = 1 / Vs), and is 
equal to 3.2 A. However, Bihamet al. (2001) found a value of 
s = 5 ■ 10 13 cm" 2 (d ? = 14 A ) for an amorphous carbon surface 
and IPerets & Biham (2006) assumed a value of s = 5 ■ 10 
cm" 3 (d s = 1.4 A) in their model of porous grains. We therefore 
treat d s as a free parameter, with values equal to 1 .4, 4.2 and 7 A. 



3.8. Multi-parameter approach 

To summarize, the GRAINOBLE model depends on eight key 
parameters presented in the previous paragraphs and listed in 
Tab. [3] along with the range and values considered in this study 
for each of them. The values reported in bold are used in the 
"reference" model (® . They correspond to the average physi- 
cal conditions of prestellar cores (density, temperature, and grain 
size), to the average values found in the litterature (site size, 
diffusion energy to binding energy ratio, initial abundance of 
atomic oxygen), to the values of parameters used by most of 
previous astrochemical models (porosity factor) or to the most 
recent values measured by experiments (activation energy). 

In total, we ran a grid of 17496 models varying these eight 
free parameters. 

3.9. Computational aspects 

The code solves three sets of differential equations: the first one 
describes the evolution of the densities of species (other than H, 
whose density is given by eq. [2J in the gas phase, the other two 
describe the chemical composition of the mantle outermost layer 
on the non-porous surface and within the pores. The equations 
giving the evolution of the chemical composition of the mantle 
layer do not depend on the grain size and are expressed in mono- 
layers/sec (MLs/s) where a monolayer is given by the number of 
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Table 3: List of the key free parameters of GRAINOBLE and the 
value range. 



Parameter 


Values 


Density n H 


10 4 - 10 s - 10 6 cm J 


Temperature T g = T d 


10 - 15 - 20 K 


Initial oxygen abundance X(0),„, 


2 • 10~ 5 - 6 • 10 5 - 2 • 10~ 4 


Grain size a c i 


0.1 -0.2-0.3yum 


Energy ratio E d /E b 


0.5 - 0.65 - 0.8 


Porosity factor F por 


- 0.3 - 0.6 - 0.9 


Site size d s 


1.4-4.2-7 A 


Activation energy E a 


400 - 1450 - 2500 K 



Notes. Bold values mark the values adopted in the reference model (see 
text, ®. 



particles of the considered species divided by the number of sites 
of a layer. The equations are the following: 



drig(i) 
dt 



= - S(i) ■ v(0 ■ n g (i) ■ cr(a d ) ■ n d {a d ) 
+ R ev (i) ■ P np (i) ■ N s (a d ) ■ n d (a d ) 



(3) 



dP np (i) 
dt 



- -S(i) ■ d s ■ v(i) ■ rig(i) - R ev (i) ■ P np {i) 

- — ■ RhopW ■ PnpW + — RhopiU ■ PporW 

** np r p 0r 



Rr,ext(if) ■ PnpUrl) ' P npihl) 



2j R r,ex,{id) ■ Pnp(irl) ' Pnpihl) 



(4) 



dPp 0r (i) 
dt 



+ ~~~ ■ Rhop(i) ■ P np (i) - -rr— • Rho P {i) ■ P por (i) 

r np ^ par 

+ ^ Rr,in(if) ' Pporihl) ' P ' porihl) 

</ 

~ R r,in(id) ■ Ppor(h l) ' P poAhl) (5) 



where P np (i) and P por (i) are the surface populations of the 
species i (in MLs) on the non-porous surface and within the 
pores respectively, and n g (i) its gas-phase density. The first and 
second terms in equation [3] describe the accretion rate and the 
evaporation rate. The first and second terms of equation [4] de- 
scribe the accretion and the evaporation rates. The third and 
fourth terms describe the exchange rate between the non-porous 
surface and the pores via diffusion. The fifth and sixth terms de- 
scribe the production and the destruction rates of i via reactions 
on the non-porous surface. Similarly, the first and second terms 
of equation|5]describe the exchange rate between the non-porous 
surface and the pores via diffusion. The third and fourth terms 
describe the production and the destruction rates of i via reac- 
tions in the pores. Note that the evolution of composition within 
all the pores can be treated with a single equation because all 
pores are assumed to have the same size. The equations giving 
the evolution of each pore can be deduced from equation [5] by 
multiplying it with a simple factor of proportionality. 



3.10. What the model does not consider 

We aim to study the impact of the multilayer behaviour and 
the porous structure of the grains on the mantle formation and 
its chemical composition. To this end, we considered thermal 
processes only. Therefore we did not consider t he desorption 
proce sses caused by exothermic surface r eactions ( Garrod et al J 
2007I). and the photolvtic desorption dGberg et al] l2009bllat 
Arasa et aT1l2010l) in spite of their relative possible importance. 
Neither did we account for the photodissociation processes on 
the ice in this version of GRAINOBLE owing to the reasons 
given in Sect. 2.1. 



4. Results 

We ran a grid of 17496 models, varying the eight free parameters 
listed in Table [3] In this section, we describe the results of our 
model. We start by describing the results relative to the reference 
model, and then we describe how the results change when each 
of the parameters changes. The analysis is based on the distribu- 
tion of the obtained mantle abundances when all parameters are 
varied, except for the parameter under consideration. This allows 
us to understand if and by how much the considered parameter 
influences the results. The plot of the multilayer and bulk distri- 
butions is shown in section l4~2l while the distributions of other 
parameters (using the multilayer approach only) are shown in 
the appendix. Finally, Table [4] summarises the main effects on 
the resulting mantle abundances caused by each parameter. 



4. 1 . The reference model 

This section describes the results obtained using the reference set 
of parameters (the boldface values in Table [3] )■ Below we refer 
to it as the "reference model". Figure [4] presents the evolution 
with time of the abundance of each species in the gas phase and 
in the grain mantles. 

The evolution of the chemical composition of the mantle 
depends on the initial gas phase abundances. In the reference 
model, O and CO abundances are approximately five times 
higher than the H abundance. Consequently, H atoms cannot hy- 
drogenate all O and CO molecules that accrete onto the mantles. 
Besides, because the CO hydrogenation reaction has a barrier of 
400 K, H atoms react preferentially with O and its products via 
barrierless reactions, ending up in water molecules. In 2 ■ 10 s 
years, 80 % of the CO reservoir is trapped within the mantle 
bulk, whereas ~ 60 % of O atoms are contained in water. The re- 
maining O atoms are shared between O3 (ozone), OH and H2O2 
(hydrogen peroxyde). Indeed, the high O to H initial gas phase 
abundance ratio (~ 6) and the relatively high temperature (15 K) 
allow a significant amout of O atoms to diffuse on the surface 
and to react with each other before meeting hydrogen atoms. 

However, gas phase CO and O abundance ratios relative to 
H gradually decrease with time (because they freeze-out onto 
the mantles), increasing the possibility of hydrogenation reac- 
tions. After ~ 5 ■ 10 4 years, the formation rates of formalde- 
hyde and methanol start to increase while the formation of ozone 
stops, corresponding to the time when H, O, and CO have sim- 
ilar abundances. The percentage of ozone, formaldehyde, and 
methanol in the ices therefore strongly depends on the age of the 
core. Young cores would show high abundances of ozone and 
of hydrogen peroxyde with respect to water, whereas older ones 
would show high abundances of formaldehyde and methanol. 

Another important result of the model is that a significant 
amount of radicals is trapped in the bulk because of the multi- 
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Time (years) 

Fig. 4: Abundance of species in interstellar grain mantles (solid) and in gas phase (in dashed) as a function of time for the reference 
model (boldface values in Table [3}. Left panel: the most abundant stable species. Middle panel: reactive species. Right panel: less 
abundant stable species. 



layer treatment. Radicals assumed to be the precursors of COMs 
such as OH, HCO or CH3O, reach abundances between ~ 5- 10~ 6 
and 10~ 9 . OH is more abundant because of its barrierless forma- 
tion reaction, and its formation mostly occurs at shorter times. 

The evolution of the chemical composition can also be stud- 
ied "spatially". Indeed, our multilayer approach allows us to 
study the composition of each monolayer within the grain man- 
tle. Figure [5] shows the evolution of the formation time of each 
layer and the evolution of the mantle thickness with time. In 
the reference model, a grain mantle of 77 layers is formed in 
2-3 ■ 10 5 yr. Assuming that the layer thickness is equal to the site 
size (here 4.2 A), it gives a total mantle thickness of ~ 0.04 fim, 
namely 32 % of the grain radius. The first 60 layers are created 
in less than 10 5 yr because they have a fast formation rate (each 
layer is formed in less than ~ 2 ■ 10 3 yr). The formation time of 
layers increases sharply at ~ 10 5 yr, because of the drop of the 
gas phase abundances. 




20 40 60 

ML number from refractory grain surface 

Fig. 5: Formation time of each monolayer (blue line) and mantle 
formation time versus its thickness expressed in monolayers (red 
line) for the reference model (boldface values in Table[3]). 



The differentiation of the composition between the inner and 
the outer layers of the mantle can also be seen in Fig. [6] Indeed, 
"intermediate" molecules such as O3, which is created from 
atomic oxygen, are rapidly formed and are therefore abundant 
in the inner layers of the mantle (the first 50 layers), whereas 
formaldehyde and methanol are mostly formed later, their frac- 
tional composition showing significant abundances only in the 
20 outermost layers. The abundances of the two main species, 
CO and H2O, remain relatively constant throughout the mantle 
up to the last outermost layers, where they decrease in favour 
of H2CO and CH3OH. The composition of radicals within the 
mantle follows the same evolution as their precursor molecules: 
the HCO abundance stays relatively constant until the drop in the 
last layers like CO, while the OH abundance decreases follow- 
ing the O abundance behaviour. Finally, the CH3O abundance 
increases approaching the surface, as H2CO. 



4.2. Multilayer versus bulk approach 

The introduction of the multilayer approach modifies the chem- 
ical behaviour of grain mantles and thus their overall chemical 
composition. Figure|7]shows the evolution of the mantle compo- 
sition with time for the reference model, adopting a multilayer 
(described in sec. 13.31 ) or a bulk (old method described in 2.1) 
approach, respectively. The figure shows that at short timescales 
(< 5 x 10 4 yr), both approaches give similar abundances for the 
main species because the H abundance is initially low compared 
to the CO and O abundances ( 94. 11 1. In contrast, O and its as- 
sociated "intermediate" molecules (OH, O2, and HO2) are more 
efficiently burned in O3, H2O2, and H2O in the bulk compared 
to the multilayer approach. 

At later times (> 5 x 10 4 yr), reactive species, which are 
trapped by the multilayer approach, continue to react in the 
bulk method. Consequently, in the bulk approach, radicals (OH, 
HCO, HO2, CH3O) and "intermediate" stable species (O, O3, 
H2O2, O2) are totally burned in less than ~ 10 6 yr. In contrast, 
the abundances of these species do not evolve after ~ 3 ■ 10 5 
yr with the multilayer approach. Finally, the predicted methanol 
and formaldehyde abundances substantially diverge in the two 
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ML number from refractory grain surface 
Fig. 6: Fractional composition of each mantle monolayer for the reference model (boldface values in Table[3]). 



approaches at times longer than ~ 10 5 yr, namely the presumed 
ages of the prestellar cores. The bulk method predicts a final (at 
10 6 yr) methanol abundance higher by factor 6 compared to the 
multilayer approach, and a formaldehyde abundance lower by a 
factor higher than 10 4 . 

To understand how robust these results are, we considered 
the whole grid of models described in 33.81 and built the distri- 
bution of the predicted abundances. Figure[8]shows this distribu- 
tion for key stable species (CO, H 2 0, H 2 CO, and CH 3 OH) and 
the radicals (OH, HCO, and CH3O) for the bulk and multilayer 
approaches for two times (10 5 and 10 6 yrs). The plots of Fig. 
[8] tell us that at short times (10 5 yr) the predicted species abun- 
dance distributions are quite similar in the bulk and multilayer 
approaches, with the exceptions of H2CO and CH3O, for which 
the distribution are slightly different. In contrast, at long times 
(10 6 yr), the distribution of all radicals are substantially different 
in the two approaches, demonstrating that the age is a key param- 
eter. The plots also provide another important information. For 
example, considering the multilayer approach results, CO and 
OH have quite peaked distributions and their predicted abun- 
dances depend slightly if at all, on the assumed values of the 
parameters: they are robust predictions. On the contrary, H2CO, 
CH3OH, HCO and CH3O have broad abundance distributions 
implying that their predicted abundances are very sensitive to 
the values of the other model parameters. 



4.3. Porous versus non-porous grains 

The influence of the grain porosity ( ^3.2t is shown in Figure [7] 
for the reference set of parameters. The presence of the poros- 
ity in the grain has the effect of slightly speeding up and in- 
creasing (by less than a factor 4) the formation of formalde- 
hyde, methanol, and CH3O, as well as to increase the destruc- 
tion of HO2, O2 and O. The reason is that porous grains trap 
H atoms more efficiently, increasing their number on the grains, 
but in the end, the H atom numbers is still limited by their low 
gas phase abundance (Fig. |7J. Furthermore, the porosity only 
slightly modifies the distribution of the results shown in Fig. 
A.l.a. Therefore, the porosity only enhances the formation of 
hydrogenated molecules by a few factors at most, regardless of 
the values of other parameters. 
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Fig. 8: Distribution of the predicted mantle abundance X mant i e of 
the key species at 10 5 (left panels) and 10 6 yr (right panels). 
The thin blue and thick red lines refer to the bulk and multilayer 
approach. The distribution has been built by considering all the 
~ 1 8000 runs of the grid (JET 
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Fig. 7: Abundance of species in interstellar grain mantles as a function of time for the reference model (solid lines), the bulk approach 
model ( {34.2I dashed lines), and the porous grain model ( {34.31 dotted lines). Left panel: the most abundant stable species. Middle 
panel: reactive species. Right panel: less abundant stable species. 



4.4. Influence of physical conditions 

Density 

The density plays an important role for two reasons. First, 
increasing the density reduces the time needed to form each 
monolayer (because the accretion rate is proportional to the 
square of the density). Second, the ratio between the number 
of H atoms and the number of heavy particles (CO and O) that 
land on the grain is inversely proportional to the density (see 
section l3~5l l. As a consequence, the higher the density, the lower 
the abundance of stable species (H2O, H2CO, CH3OH) created 
by the hydrogenation reactions, especially at long times (> 10 s 
yr) (Fig. A.l.b). In contrast and for the same reasons, mantle 
CO and OH are more abundant for higher density. 

Temperature T g — Tj 

The grain temperature has a moderate influence on the 
predicted mantle composition, although the diffusion (and thus 
the reaction) and the desorption rates depend exponentially on 
it. However, the desorption rate increases much more quickly 
than the diffusion rate with temperature (by factors of 10 2 - 10 6 
depending on the E^/Ei, ratio). Consequently, with increasing 
grain temperature, H atoms have a higher probability to desorb 
into the gas phase before encountering another particle. The 
final mantle abundance of molecules created by hydrogenation 
reactions (H2CO and CH3OH) is therefore lower at any time 
at higher temperatures, as shown by Fig. A. I.e. The mantle 
abundance of radicals, on the other hand, is not affected by the 
grain temperature. 

Grain size aj 

While the mantle thickness depends roughly quadratically 
on the grain size, the final mantle composition depends on it 
very little. The reason is that the particle accretion rate (inversely 
proportional to the grain size square) and the formation time of 
each layer compensate each other. The result is that the mantle 
composition does not depend on the assumed grain size, whether 
0.1 or 0.3 pm. 



4.5. Influence of other key parameters of GRAINOBLE 

Initial abundance of atomic oxygen X( 0) ml 

Not surprisingly, the initial gas phase abundance of atomic 
oxygen is an important parameter for the final mantle water 
abundance. The higher the O abundance, the higher the iced 
H2O. Less evidently, oxygenation is in competition with hydro- 
genation. The initial O abundance also affects the final mantle 
abundance of formaldehyde and, even more, methanol, because 
the elemental C/O ratio decreases between 0.7 and 0.2 when 
X(0)i„i increases between 2 ■ 10~ 5 and 2 • 10~ 4 . At 10 s yr, the 
predicted methanol abundance is a factor ten higher when the O 
abundance is a factor ten lower. 

Diffusion energy to binding energy ratio E c i/Eb 

The diffusion energy is an important parameter in the 
predicted mantle abundance of formaldehyde and methanol. 
Because the diffusion rates depend exponentially on the dif- 
fusion energy, the mantle abundance of stable species, like 
formaldehyde and methanol, decreases with increasing E^. 
Conversely, the diffusion energy does not affect the mantle 
abundance of radicals, because their formation and destruction 
rates compensate each other. 

Activation energy E a 

The activation energy of CO and H2CO hydrogenation re- 
actions strongly influences the formation rates of formaldehyde 
and methanol because their reaction rates depend exponentially 
on E a . Most of the runs with a low value of E a (400 K) predict 
mantle abundances of formaldehyde and methanol higher than 
10~ 6 at > 10 5 yr. Conversely, higher values of E„ give uniform 
distributions of H2CO and CH3OH abundances, between 10~ 18 
and 10~ 6 . Only ~ 5% of the runs using a activation energy equal 
to 2500 K predict formaldehyde and methanol abundances 
higher than 10~ 6 . The final abundances of formaldehyde and 
methanol, consequently strongly depend on the other model 
parameters for high E a values. A significant amount of radicals, 
like HCO and CH3O, can survive in the mantle only if the 
activation energies are low. Runs with activation energies of 
400 K predict radical mantle abundances between 10" 10 and 
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Table 4: Summary of the effect of each free parameter on the surface chemistry and on the distribution of results (shown in the 
appendix). 



Parameter 


Effect on the surface chemistry 


Distributions of abundances for stable 


Distributions of abundances 






species 


for radicals 


Multilayer 


The multilayer approach 


- X(CO) > 10- 5 for 


- X(OH) > 10-" for 


versus bulk 


- decreases the formation of stable 


88% (ML), 50% (bulk) cases 


75% (ML), 1% (bulk) cases 


approach 


species 


- X(H 2 0) > 1(T 5 for 


- X(HCO) > 10- 9 for 




- traps reactive species in the mantle 


60% (ML), 93% (bulk) cases 


25% (ML), 2% (bulk) cases 






- X(CH 3 OH) > 10~ 5 for 








15% (ML), 55% (bulk) cases 




Porosity 


The increase of porosity 


-X(H 2 CO) > 10-' for 


- X(OH) > 10-" for 




- traps volatile and diffusing particles 


52% (smooth), 61% (porous) cases 


55% (smooth), 20% (porous) cases 




- increases the formation rate of stable 


- X(CH 3 OH) > 10 - 1 for 






species 


51% (smooth), 61% (porous) cases 




Density 


The increase of density 


- X(H 2 0) > 10" 3 for 


-X(OH) > -lO" 6 for 




- decreases the formation of stable 


2% (dense) , 85% (sparse) cases with 


0% (sparse), 75% (dense) cases 




species formed by hydrogenation 


X(H 2 0) -> 2 • 10~ 4 






reactions 


- X(CH 3 OH) > 10~ 5 for 








0% (dense) , 35% (sparse) cases with 








X(CH 3 OH) -» 9 • 10- 5 




Temperature 


The increase of temperature 


-X(CH 3 OH) > 10"' for 


- X(OH) > 10~ b for 




- strongly increases the desorption of H 


65% (cold), 45% (warm) cases 


55% (cold), 48% (warm) cases 




andO 


- X(H 2 CO) > 10~ 7 for 






- decreases the rates of hydrogenation 


60% (cold), 45% (warm) cases 






reactions 








- decreases the abundance of stable 








species formed by these reactions 






Grain size 


The increase of the grain size 


Same evolution of distributions for all 


Same evolution of distributions for all 


CLd 


- increases the mantle thickness 


stable species 


radicals 




- does not influence the integrated 








mantle composition 






Initial 


The increase of X ini (0): 


- X(H 2 0) > 2 • 10- 5 for 


- X(OH) > 4 • lO 6 for 


abundance of 


- increases the formation of water (and 


0% (low-O), 65% (high-O) cases with 


0% (low-O), 50% (high-O) cases with 


atomic oxygen 


other molecules formed from reactions 


X(H 2 0) -> 2 ■ 10- 4 


X(OH) -> 4 • 10~ 5 




involving O) 


- No influence is seen for CO, H 2 CO 






- slightly decreases the formation of 


and CH 3 OH 






molecules formed from CO 






Diffusion 


The increase of the diffusion energy 


- X(H 2 0) > 5 • 10- 6 for 


Same evolution of distributions for all 


energy 


- decreases the diffusion rate of mobile 


87% (fast), 74% (slow) cases 


radicals 




species 


- X(CH 3 OH) > 10~ 7 for 






- decreases the formation rate of stable 


80% (fast), 30% (slow) cases 






species 


- X(H 2 CO) > 10~ 7 for 






- slightly increases the survival of 


75% (fast), 40% (slow) cases 






radicals 






Activation 


The increase of activation energies 


- X(CO) > 10 5 for 


- X(HCO) > 10 1U for 


energies 


- strongly decreases the formation rate 


35% (low-Ea), 98% (high-Ea) cases 


90% (low-Ea), 0% (high-Ea) cases 




of H 2 COandCH 3 OH 


- X(H 2 CO) > 10- 7 for 


- X(CH 3 0) > 10 10 for 




- strongly decreases the survival of 


83% (low-Ea), 25% (high-Ea) cases 


70% (low-Ea), 0% (high-Ea) cases 




radicals HCO and CH 3 


- X(CH 3 OH) > 10~ 6 for 






- slightly increases the reaction rates of 


80% (low-Ea), 5% (high-Ea) cases 






reactions involving O 








- slightly increases the formation of 








water 






Site size 


The increase of the site sizes 


Same evolution of distributions for all 


Same evolution of distributions for all 




- decreases the formation time of layers 


stable species 


radicals 



- increases the mantle thickness 

- does not influence the integrated 
mantle composition 



10 for both radicals, while runs with E a = 2500 K predict 
abundances less than ~ 10~ n . Not surprisingly, the activation 
energy of CO and H 2 CO hydrogenation reactions are, there- 
fore, critical parameters for the survival of radicals in the mantle. 

Site size d s 

The reaction rates and the formation time of each monolayer 
are both functions of the site size d s : the former is proportional 



to d 2 s , whereas the latter is proportional to dj 2 . Therefore, the 
two processes, which are in competition, cancel each other, and 
consequently, the total mantle abundance of stable species and 
radicals do not sensitively depend on the size of sites. 

However, the thickness of the grain mantle strongly depends 
on the site size. Indeed, the maximum number of monolayers is 
given by the ratio between the number of particles on the grain 
N par t and the number of grain sites N s (proportional to d~ 2 ). 
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Furthermore, if we assume that the thickness of a layer is equal 
to the distance between two sites, the thickness of the mantle is 
proportional to dj. Thus, the bigger are the sites, the faster the 
grains grow. Consequently, the depletion is more efficiently on 
grains with big sites, increasing the abundance of methanol and 
decreasing the abundance of CO for example. 

4.6. Concluding remarks 

Based on the previous paragraphs, we can identify three classes 
of parameters, depending on their influence on the results: 

- 1) The parameters that have a strong influence on the man- 
tle composition are the approach to model the chemical be- 
haviour (multilayer versus bulk), the density, the diffusion 
energy and the activation energy of CO and H2CO hydro- 
genation reactions. Modifying the values of these parameters 
within their chosen range drastically changes the distribution 
of the predicted mantle abundances, both in the shape and in 
the values. 

- 2) The grain porosity and the temperature as well as the ini- 
tial abundance of atomic oxygen have only a moderate im- 
pact on the predicted mantle composition. 

- 3) The grain and the site sizes have a negligible influence 
on the total mantle composition. However, these parameters 
strongly modify the ice thickness, which could have an im- 
pact on photolytic processes and for the desorption time of 
the mantle during the warm-up phase in protostar envelopes. 

5. Comparisons with previous microscopic models 

The Monte-Carlo continuous ti me random w alk (C TRW) 
method was i nt roduce d by IChang et al.l (|2005) and 
Cupp en & Herbsll d2005l) for the hydrogen recombination 
system, and was th e n extended to mor e comp le x networks in 
Chan g et al.l (120071) . ICuppen & Herbstl (120071) . ICuppen et all 
(I2009h . The natural output of this model is the multilayer 
structure of the mantle. This section compares our results 
with the model of Cup pen et alJ d2009h with the twofold aim 
of validating the GRAINOBLE code and of highlighting the 
differences between the two models. 

5.1. Validation of the GRAINOBLE code 

Cup pen et al.l d2009l) considered the accreti on of H, CO, and 
H2 on smooth surfaces as described by ICuppen & H erbst 
(2005), leading to the formation of formaldehyde, methanol, 
and their associated radicals HCO and CH3O via the Langmuir- 
Hinshelwood and the Eley-Rideal mechanisms. A major differ- 
ence between GRAINOBLE and the Cuppen et al. model is that 
they consider an a-CO ice, whereas we assumed that the bulk of 
the mantle is formed by iced water. To reproduce the interactions 
of molecules within an a-CO ice, we modified the set of binding 
energies in our model accordingly. Second, the CTRW model 
takes into account the individual interaction energies between 
each particle, while our model only considers the total binding 
energy between the adsorbate and the substrate. To simulate their 
approach, we assumed for each species i a binding energy equal 
to fourteen times the interaction between the species i and a CO 
molecule (£/,(/) = WZsy-coX corresponding to the interaction 
between an adsorbate and a porous multi-layered CO ice. The 
binding energy of H is equal to 450 K, while the binding energy 
of CO is 880 K. The di ffusion to desorption energy ratio is as- 
sumed to be 0.8, as in ICuppen et al. (2009) for a layer that is 



initially smooth. The activation energies fo r CO and H2CO hy - 
drogenation reactions are those measured bv lFuchs et alJ {2009), 
and slightly depend on the grain temperature. 

To validate ou r code with respect to the model of 
Cup pen et alJ (I2009I) . we run our model with a density of 10 5 
cm -3 and considered the res ults at 2 X 10 s yr, fo r the following 
four cases (as in Part 3.3 of ICuppen et al.l d2009)): a) grain tem- 
perature 12 K, n(H)/n(CO)=l and n d = 2 x \0 n n h ; b) as a) 
but n L t — 1 x 10 12 n/,; c) as a) but temperature 15 K; d) as b) but 
n(H)/n{CO) =0.5. Figure [9] shows the chemical composition of 
each monolayer for the four cases, for porous grains (F por = 0.7: 
this is the value that agrees best for o ur and Cuppen's mo del). 
The comparisons between Figure 4 of Cup pen et al.l d2009t) and 
our Figure [9] shows that the two codes give very similar results, 
which validates our code. 




ML number from refractory grain surface 

Fig. 9: Chemical composition of the mantle layers as a function 
of the mon olayers at a time of 2 - 10 s years, for a density h# = 10 5 
cm 3 , as in Cuppen et al.l d2009t) : see text): CO (red), formalde- 
hyde (blue) and meth anol (green). The fou r panels refer to the 
four cases of Fig.4 of [Cuppen et alJ d2009) (see text): a) grain 
temperature 12 K, n{H)ln{CO)-\ and n d — 2 x 10 12 «/, ; b) as a) 
but rid — 1 x 10 12 n/,; c) as a) but temperature 15 K; d) as b) but 
n(H)/n(CO)=0.5. 



5.2. Differences between the two approaches 

We highlight a few critical differences between the model of 
Cuppen and ours. 

Unlike our multilayer approach, the CTRW method allows 
us to model the segregation effects. Heavier particles can ag- 
glomerate together to form islands where they can be fixed 
onto the grains more efficiently. Our multilayer method is a 
macroscpic approach where segration effects are not taken into 
account. However, the pores can play the same role as the is- 
lands, because particles that enter into the pores cannot desorb 
either. This is supported by the fact that our porous grains case 
agrees best with Cuppen's model. 

An important difference between ours and Cuppen et al.'s 
model is the inclusion of the water ice formation. Because 
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Cuppen et al. did not consider the formation of water, their 
ices are formed by pure a-CO ice. Consequently, the interaction 
energies are lower than those used in our model. Therefore in 
Cuppen's model ices can form only at temperatures lower than 
16-18 K, whereas in our model iced CO can still survive un- 
til 21-23 K. In our model, then, the formation of formaldehyde 
and methanol can oc cur in a wider range of temperatures than in 
Cup pen et al.l 12009). Consequently, the abundance of formalde- 
hyde and methanol formed on grain mantles decreases much 
more efficiently with the increase of the temperature in Cuppen 
et al.'s model than in ours. For example, using the same physical 
parameters, the CTRW model predicts a decrease of the mantle 
thickness of 85 % (~ 40 to 6 monolayers), whereas our model 
does not predict any thickness decreases between 15 and 16.5 K. 

6. Comparisons with observations 

6.1. The observations 

In order to compare our model predictions with observations, we 
used the data pertaining to the ices only, because the abundance 
of the sublimated mantle species does not necessarly reflect the 
chemical composition of their precursor ices. First, it can be al- 
tered by reactions in the gas phase. Second, the different subli- 
mation temperatures may introduce errors when comparing the 
abundance ratios of different species, because they may refer to 
different regions. 

Unfortunately, solid CO, H 2 CO, and CH 3 OH have been 
observed simultaneously only towards a v ery small sample 
of high- and intermed iate-mass protostars dGibb et al.l [2004; 
IPonto ppidan et al. 2004). To increase the statistics and to include 
low-mass protostars, too, we therefore restricted the comparison 
of our model to the observed solid CO and CH3OH abundances 
only. Fortunately, a posteriori, when considering the model pre- 
dictions, we found that using the H2CO would not provide (sub- 
stantially) more constraints. Table [5] reports the compilation of 
observations that we used for the comparison. 



Table 5: Mantle abundances of CO and CH3OH with respect to 
water along with the relative CH3OH/CO abundance ratio, as 
observed towards low- (LM), intermediate- (IM), and high-mass 
(HM) protostars. 



Source 


X(CO) 


X(CH 3 OH) 


X(CH 3 OH) 




% wrt H 2 


% wrt H 2 


% wrt to CO 


Low-Mass 








LI 448 IRS 1 


45.5 a -" 


< 14.9 d 


< 32.7 


L1455 SMM 1 


5.6 a - e 


< 13.5 d 


< 241.1 


RNO 15 


13.6 *-' 


< 5.0 d 


< 36.8 


IRAS 03254 


13.6 °- e 


<4.6 d 


< 33.8 


IRAS 03271 


8.2 


< 5.6 d 


<68.3 


Bl-a 




< 1.9 d 


< 14.6 


Bl-c 


28.6 "- e 


< 7.1 d 


< 24.8 


L1489 


15.2 °- e 


4.9 d 


32.3 


DG Tau B 


20.5 


< 5.7 d 


< 27.8 


IRAS 12553 


12.6 


< 3.0 d 


<23.8 


IRAS 13546 


22.7 "~ e 


< 3.9 d 


< 17.2 


IRAS 15398 


40.3 a -" 


10.3 d 


25.6 


CRBR 2422.8-3423 


11.2 


< 9.3 d 


< 83.0 


RNO 91 


19.0 "- e 


< 5.6 d 


< 29.5 


IRAS 23238 


4.0 


< 3.6 d 


90.0 


Intermediate-Mass 








AFGL989 


19.0 * 


1.7 * 


8.9 


SMM4 


30.0 c 


28.0 c 


93.3 


High-Mass 








W33A 


7.4 "- e 


14.7 d 


198.6 


GL 2136 


10.2 "-" 


8.5 d 


83.3 


S140 


16.6 "- e 


< 3.0 d 


< 18.1 


NGC 7538 IRS 9 


16.5 


7.5 d 


45.5 


AFGL2136 


5.0 6 


5.0* 


100.0 


AFGL7009 


16.0 6 


33.0* 


206.3 



Refe rences. (a) IPontoppidan et alj d2003l); m iGibbet alJ ([2 004); 
<c) IPontoppidan et al.1 d2004h : w iBooeert et alJ d20081)~ < e) 
IPontoppidan et al.l ( 120081) . 



6.2. Constraining the activation and diffusion energies 

As a first step, we compared the model predictions with the ob- 
servations with the goal to constrain the two microphysics pa- 
rameters, which have a high impact on the model predictions 
(§4.6), namely E d /E b and E a . Figure [TOl shows the CH3OH/CO 
abundance ratio as a function of time obtained assuming refer- 
ence parameters (and more particularly a density 10 5 cm -3 ) and 
for different values of the activation energy E a and diffusion en- 
ergy E4. The comparison of the model predictions with the val- 
ues observed towards intermediate- and high- mass protostars 
clearly excludes an activation energy higher than 1450 K and a 
diffusion energy higher than 0.65 times the binding energy E/,. 
More stringently, if E a = 1450 K, Ej/Eb has to be equal to 0.5. 
The analysis of a similar plot for a density of 10 4 cirT 3 gives 
similar constraints: £0=1450 K implies Ed/Eb-0.5 and £0=400 
K, Ed/Et, ~ 0.65. The case 10 6 cm -3 does not provide more con- 
straints, it is indeed excluded by the observations (see below). 



6.3. Constraining the pre-collapse phase duration and 
density 

Once we limited the range of possible values of E a and E\, con- 
sistent with the observations, we can attempt to constrain the 
duration of the pre-collapse phase and the density. To this end, 



Figure QT| shows the curves of the CH3OH/CO abundance ra- 
tio versus time for three different densities (10 4 , 10 s and 10 6 
cm -3 ) and computed for the two following sets of E a and Ej 
values: 1) £,,=400 K and Ed/Eb=0.65 (and the reference param- 
eters) ; 2) £0=1450 K and Ed/Eb-Q.5 (and the reference pa- 
rameters). The comparison between the predicted and observed 
CH3OH/CO abundance ratio in intermediate- and high-mass 
protostars suggests that the bulk of methanol has been formed 
when the pre-collapse condensation had densities between 10 4 
and 10 5 cm -3 . Higher densities would produce an insufficiently 
CH3OH/CO abundance ratio compared to the values observed 
in intermediate- and high- mass protostars. If the densities are 
around 10 5 cm -3 , the duration of the pre-collapse phase must 
have lasted at least 10 5 yr. Lower densities would allow shorter 
pre-collapse duration times. But densities higher than 10 s cm -3 
are not ruled out for the condensations that produce the low-mass 
protostars. In that case, the duration of the pre-collapse phase 
lasted more than 5 x 10 4 yr. 

One important result is that the CH3OH/CO abundance ratio 
does not provide an estimate of the duration of the pre-collapse 
phase, but only a lower limit to it. It neither provides an esti- 
mate of the density at the time of the collapse but only a range 
of possible densities for any given observed CH3OH/CO abun- 
dance ratio. Note that the CH3OH/H2CO abundance ratio does 
not allow us to arrive at better constraints, because the theoretical 
curves are very similar to those of Fig.fTTI 



14 



V. Taquet et al.: The GRAINOBLE model 



10 £ 



E d /E b =0.5 

E d /E b =0.65 

E„/E h =0.8 



solid: E a = 400 K 
dotted: E a = 1450 K 
dashed: E a = 2500 K 
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Fig. 10: Mantle CH3OH/CO abundance ratio versus time. The 
blue, green, and red lines refer to Ed/Et,-0.65, 0.5, and 0.8, re- 
spectively. Solid, dotted, and dashed lines refer to £,,=400 K, 
1450 K, and 2500 K, respectively. In these computations the 
density is 10 5 cm 3 . The box with hatching shows the interval 
of CH3OH/CO abundance ratio observed towards intermediate- 
and high- mass protostars (Table [5J. Because low-mass proto- 
stars provide only less stringent upper limits, they are not re- 
ported in the plot. 



10' 



Model 1 solid: n H = 10* cm 3 
Model 2 dotted: n H = 1 0' cm 3 
dashed: n H = 10 ! cm 3 



o 10 1 
o 




10 4 

Time (years) 

Fig. 11: Mantle CH3OH/CO abundance ratio versus time. The 
red lines refer to Model 1 (£ fl =400 K and E d /E b =0.65), the blue 
lines refer to Model 2 (£ a =1450 K and E d [E b =Q.5) (see text). 
Solid, dotted and dashed lines refer to densities 10 4 , 10 5 and 
10 6 cm -3 respectively. The grey box with hatching shows the 
interval of CH3OH/CO abundance ratio observed towards high- 
and intermediate- mass protostars, while the green dashed box 
shows the values (upper limits) observed in low-mass protostars 
(Table©. 



6.4. Chemical differentiation within grain mantles 

As briefly mentioned in §2, IR observations suggest a chemical 
differentiation of the grai n mantles with a p olar matrix formed 
already at A v ~ 3 mag dWhittet et al]ll988l). and a non- polar 
CO matrix formed at higher densities dTielens et al.lll99lh . In 
between, CO can b e mixed either with water, methanol, or both 
at the same time (Bisschop et al. 2007), as well as mixed with 
O2, O3, N 2 , or CO2 (Pontoppidan 2006). Actually, this observa- 
tional fact is one of the motivations of our work. As discussed 



in §|4]and, for example, shown in Fig. [6] our code predicts inho- 
mogeneous grain mantle composition. For the reference model, 
for example, CO is mixed with water in the inner layers while 
a mixture between CO and methanol is predicted in the outer 
layers. Given the lack of a dataset to compare our results with, 
it is difficult to carry out quantitative comparisons between ob- 
servations and predictions. However, it is worth noticing that the 
exact grain mantle composition and stratification depends on the 
overall evolution and structure of the pre-collapse condensation, 
because the chemical differentiation depends on the density, the 
temperature and time. 

7. Conclusions and perspectives 

We presented a new model, GRAINOBLE, which computes the 
chemical composition of interstellar grain mantles during the 
cold and dense phase of the pre-collapse. The model presents 
two main differences from most other published codes: 1) it as- 
sumes that only the outermost mantle layer is chemically active, 
while the mantle bulk is not; 2) it considers porous grains. 

We run GRAINOBLE for a large set of parameters, which 
are either unconstrained (density, temperature and grain size of 
the prestellar condensation) or poorly known (gas atomic oxy- 
gen abundance, diffusion energy, hydrogenation reactions activa- 
tion energy, grain site size) . We obtained a grid of about 18000 
models and built distribution plots of the predicted mantle abun- 
dances. This allows us to study the influence of each of these 
parameters on the predicted abundances and, consequently, the 
robustness of the predictions. 

Finally, we compared the GRAINOBLE predictions 
with those obtaine d by the microscopic CTRW model of 
Cupp en et al.1 (120091) . obtaining a fair agreement between the 
two models. The comparison clearly validates our treatment and 
our code (§0. The most important results of the GRAINOBLE 
model are the following 

1) The multilayer treatment shows a differentiation of the species 
within the mantle, with the innermost layers being rich in CO, 
the intermediate layers in formaldehyde, and the outermost lay- 
ers in methanol, reflecting the different formation time of each 
species in the mantle. This differentiation will likely lead to a 
differentiation in the deuteration of formaldehyde and methanol 
on the ices. Because t he deuteration increases with the CO deple- 
tion i n the gas phase dRoberts et al.ll2003l:ICeccarelli & Dom inik 
2005), methanol will likely be more deuterated than formalde- 
hyde, as indeed observed in lParise et alj d200 6). A forthcoming 
paper will focus on the detailed modeling of the deuteration pro- 
cess. 

2) The multilayer treatment predicts a relatively high abundance 
of radicals trapped in the mantle. For example, HCO and CH3O 
show abundances of 10 - 10 with respect to H nuclei. As 
suggested by Garrod & Herbst (2006), these radicals can react 
to form COMs when the grain temperature increases during the 
formation of the protostar. At this stage, we cannot say whether 
the predicted radical abundances are sufficient to explain the 
present observations (e.g. Fig. 1), and, unfortunately, compar- 
isons with previous models are not possible because of the lack 
of specific information. We will explore this aspect, the forma- 
tion of COMs, in a forthcoming paper. 

3) The presence of porosity in the grains only moderately in- 
fluences the mantle chemical composition, causing an enhanced 
abundance of formaldehyde and methanol by less than a factor 
four. 

4) The chemical composition of grain mantles strongly depends 
on the physical conditions of the prestellar condensation, par- 
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ticularly on its density as well as its age. Therefore, it will be 
important to model the evolution of the condensation for real- 
istic predictions of the mantle composition. Conversely, the ob- 
served mantle composition can provide valuable information on 
the past history of protostars, and will be the focus of a forthcom- 
ing paper. Comparisons of the present predictions with the ices 
observations (specifically the CH3OH/CO abundance ratio) sug- 
gest that intermediate- and high- mass protostars evolved from 
condensations less dense than about 10 5 cm 4 , while no strin- 
gent constraints can be put on the prestellar condensations of 
low-mass protostars. 

5) The predicted mantle composition critically depends on the 
values of the diffusion energy and activation energy of the hy- 
drogenation reactions. Comparing our model predictions with 
observations of ices allows to constrain the values of these two 
important parameters: the diffusion to binding energy ratio has 
to be around 0.5-0.65, and the activation energy has to be less 
than 1450 K. 

6) Other parameters of the GRAINOBLE model, like the gas 
atomic oxygen abundance and the site size, do not substantially 
influence the predicted mantle abundances. 

In summary, GRAINOBLE is a versatile and fast code to 
model the grain surface chemistry. Its main improvement with 
respect to previous similar models is the treatment of the layer- 
by-layer chemistry, which provides a more realistic stratified 
mantle composition. It predicts that the mantles are indeed strat- 
ified and that radicals are trapped inside the mantles. Because it 
is not expensive from the point of view of computing time, mod- 
eling of complex and more realistic cases, like the evolution of a 
prestellar core, are now feasible. 
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Appendix A: Distributions of abundances 

As introduced in section 3.8, we have built a model grid in which 
eight input parameters vary. To study the results of this grid, 
we computed the distributions of selected species abundances on 
grain mantles X m (i), namely the number of results giving X m (i) 
in each abundance interval. The following figures show the dis- 
tributions of abundances for each input parameter and for two 
extremal values, using the multilayer approach only. 
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Fig. A.l: Distribution of abundances on grain mantles at t — 10 5 
and 10 6 years for two porosity values: F por = (thick red) and 
0.9 (narrow blue). 



Appendix B: List of symbols used in this work 
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Table B.l: List of symbols and their explanation. 



Symbol 


Parameter 


ad 


diameter of interstellar grains (in /im) 


d s 


size of each site (in A ) 


E d 


barrier energy against diffusion of a considered species (in K) 


Ei, 


barrier energy against desorption of a considered species (in K) 


E a 


activation energy of the CO and H 2 CO hydrogenation reactions (in K) 


F 

1 por 


fraction of the grain surface occupied by the pores 


F 

1 up 


fraction of the grain surface occupied by the non-porous surface 


Fed 


fraction of the grain surface occupied by the edge porous sites 


?d 


radius of interstellar grains (in /mi) 


n d 


density of interstellar grains (in cirT 3 ) 


N 

1 ' pore 


number of sites occupied by each pore 


N s 


number of sites of the considered layer 


I'M 


total density of H nuclei (in cirT 3 ) 


n i 


gas phase density of a considered species (in enr 3 ) 


P 

1 np 


surface population of a considered species in the non-porous surface (in monolayers) 


P 

1 por 


surface population of a considered species in the pores (in monolayers) 


R-acc 


accretion rate of a considered species (in s _1 ) 


Rfmp 


hopping rate from site to site of a considered species (in s ) 


Rdi/f 


diffusion rate of a considered species (in s -1 ) 




reaction rate of a considered reaction in the non-porous surface (in s ) 


Rrpor 


reaction rate of a considered reaction in the pores (in s -1 ) 


R er 


evaporation rate of a considered species (in s -1 ) 


S 


sticking coefficient of a considered species 


(Til 


cross section of interstellar grains (in cm 2 ) 


T 

1 s 


gas temperature (in K) 


T d 


dust temperature (in K) 


V 


thermal velocity of a considered species in the gas phase (in cm s _1 ) 


X,, 


abundance of interstellar grains relative to H nuclei 


X, 


abundance relative to H nuclei of a considered species in the gas phase 




abundance relative to H nuclei of a considered species on grain mantles 


e 


cosmic rays ionization rate (in s -1 ) 



3000 
2500 
2000 
1500 
1000 
50 g 

1500 
1000 
500 


400 
300 
200 
100 


500 
400 
300 
I 200 
i 100 
■ 

800 
600 
400 
200 




250 
200 
150 
100 

5 8 



CO 



1 B years 



H 2 CO 



CH.OH 



OH 




HCO 



IU ■ 



CH 3 




1 6 years 



-1 




10"'° 10 s 10' 6 10" 10 10 10" B 10" 6 10' 4 



2000 
1500 
1000 
500 


800 
600 
400 
200 


300 
250 
200 
150 
100 

5 8 

250 
200 
. 150 
! 100 

i 5 8 

500 
400 
300 
200 
100 

600 
500 
400 
300 
200 
101 




150 
100 
50 



CO 



1 5 years 



A 



H z CO 



CH 3 OH 





La 



ch,o 





10° years 






10" B 10" 6 10" 10"'' 
Y 



Fig. A. 2: Distribution of abundances on grain mantles at t — 10 5 
and 10 6 years for two densities: n# = 10 4 (thick red) and 10 6 
cirT 3 (narrow blue). 



Fig. A. 3: Distribution of abundances on grain mantles at t — 10 5 
and 10 6 years for two temperatures: T — 10 (thick red) and 20 K 
(narrow blue). 
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Fig. A. 4: Distribution of abundances on grain mantles at t — 10 5 
and 10 6 years for two grain sizes: aj — 0.1 (thick red) and 0.3 
fim (narrow blue). 



Fig. A. 6: Distribution of abundances on grain mantles at t = 10 5 
and 10 6 years for two energy ratios: Ej/Eb = 0.5 (thick red) and 
0.8 (narrow blue). 
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Fig. A. 5: Distribution of abundances on grain mantles at t - 
and 10 6 years for two initial oxygen abundances: X (m (0) 
10~ 5 (thick red) and 2 • 10~ 4 (narrow blue). 
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Fig. A. 7: Distribution of abundances on grain mantles at t = 10 5 
and 10 6 years for two site sizes: d s = 1.4 (thick red) and 7 A 
(narrow blue). 
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Fig. A. 8: Distribution of abundances on grain mantles at t = 10 5 
and 10 6 years for two CO and H2CO hydrogenation reactions: 
E a = 400 (thick red) and 2500 K (narrow blue). 
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